# Replication Archive for 
# Alexander Coppock, Emily Ekins and David Kirby (2018), 
# "The Long-lasting Effects of Newspaper Op-Eds on Public Opinion",
# Quarterly Journal of Political Science: Vol. 13: No. 1, pp 59-87.

rm(list = ls())

# Uncomment to install
#install.packages(c("tidyverse", "estimatr", "reshape2", "xtable"))

library(tidyverse)
library(estimatr)
library(reshape2)
library(xtable)

# Set working directory to replication archive

load("mturk_opeds_cleaned.rdata")
load("elite_opeds_cleaned.rdata")
load("dv_df.rdata")
source("opeds_source.R")


fit_01 <- lm_robust(dv_amtrak_agree_w1 ~ Z , data = filter(mturk_opeds, responded_w1, Z %in% c("control", "amtrak")))
fit_02 <- lm_robust(dv_climate_agree_w1 ~ Z, data = filter(mturk_opeds, responded_w1, Z %in% c("control", "climate")))
fit_03 <- lm_robust(dv_flat_agree_w1 ~ Z, data = filter(mturk_opeds, responded_w1, Z %in% c("control", "paul")))
fit_04 <- lm_robust(dv_vets_agree_w1 ~ Z, data = filter(mturk_opeds, responded_w1, Z %in% c("control", "veterans")))
fit_05 <- lm_robust(dv_wall_agree_w1 ~ Z, data = filter(mturk_opeds, responded_w1, Z %in% c("control", "wallstreet")))
fit_06 <- lm_robust(dv_amtrak_agree_w1 ~ Z, data = filter(elite_opeds, responded_w1, Z %in% c("control", "amtrak")))
fit_07 <- lm_robust(dv_flat_agree_w1 ~ Z, data = filter(elite_opeds, responded_w1, Z %in% c("control", "paul")))
fit_08 <- lm_robust(dv_vets_agree_w1 ~ Z, data = filter(elite_opeds, responded_w1, Z %in% c("control", "veterans")))
fit_09 <- lm_robust(dv_wall_agree_w1 ~ Z, data = filter(elite_opeds, responded_w1, Z %in% c("control", "wallstreet")))

fit_11 <- lm_robust(dv_amtrak_agree_25_w1 ~ Z , data = filter(mturk_opeds, responded_w1, Z %in% c("control", "amtrak")))
fit_12 <- lm_robust(dv_climate_agree_25_w1 ~ Z, data = filter(mturk_opeds, responded_w1, Z %in% c("control", "climate")))
fit_13 <- lm_robust(dv_flat_agree_25_w1 ~ Z, data = filter(mturk_opeds, responded_w1, Z %in% c("control", "paul")))
fit_14 <- lm_robust(dv_vets_agree_25_w1 ~ Z, data = filter(mturk_opeds, responded_w1, Z %in% c("control", "veterans")))
fit_15 <- lm_robust(dv_wall_agree_25_w1 ~ Z, data = filter(mturk_opeds, responded_w1, Z %in% c("control", "wallstreet")))
fit_16 <- lm_robust(dv_amtrak_agree_25_w1 ~ Z, data = filter(elite_opeds, responded_w1, Z %in% c("control", "amtrak")))
fit_17 <- lm_robust(dv_flat_agree_25_w1 ~ Z, data = filter(elite_opeds, responded_w1, Z %in% c("control", "paul")))
fit_18 <- lm_robust(dv_vets_agree_25_w1 ~ Z, data = filter(elite_opeds, responded_w1, Z %in% c("control", "veterans")))
fit_19 <- lm_robust(dv_wall_agree_25_w1 ~ Z, data = filter(elite_opeds, responded_w1, Z %in% c("control", "wallstreet")))

fit_21 <- lm_robust(dv_amtrak_agree_75_w1 ~ Z , data = filter(mturk_opeds, responded_w1, Z %in% c("control", "amtrak")))
fit_22 <- lm_robust(dv_climate_agree_75_w1 ~ Z, data = filter(mturk_opeds, responded_w1, Z %in% c("control", "climate")))
fit_23 <- lm_robust(dv_flat_agree_75_w1 ~ Z, data = filter(mturk_opeds, responded_w1, Z %in% c("control", "paul")))
fit_24 <- lm_robust(dv_vets_agree_75_w1 ~ Z, data = filter(mturk_opeds, responded_w1, Z %in% c("control", "veterans")))
fit_25 <- lm_robust(dv_wall_agree_75_w1 ~ Z, data = filter(mturk_opeds, responded_w1, Z %in% c("control", "wallstreet")))
fit_26 <- lm_robust(dv_amtrak_agree_75_w1 ~ Z, data = filter(elite_opeds, responded_w1, Z %in% c("control", "amtrak")))
fit_27 <- lm_robust(dv_flat_agree_75_w1 ~ Z, data = filter(elite_opeds, responded_w1, Z %in% c("control", "paul")))
fit_28 <- lm_robust(dv_vets_agree_75_w1 ~ Z, data = filter(elite_opeds, responded_w1, Z %in% c("control", "veterans")))
fit_29 <- lm_robust(dv_wall_agree_75_w1 ~ Z, data = filter(elite_opeds, responded_w1, Z %in% c("control", "wallstreet")))

fitlist <- list(fit_01, fit_02, fit_03, fit_04, fit_05, fit_06, fit_07, fit_08, fit_09,
                fit_11, fit_12, fit_13, fit_14, fit_15, fit_16, fit_17, fit_18, fit_19,
                fit_21, fit_22, fit_23, fit_24, fit_25, fit_26, fit_27, fit_28, fit_29)


fit_df <- 
  fitlist %>%
  map_df(tidy, .id = "fit") %>%
  filter(term != "(Intercept)") %>%
  mutate(study = rep(rep(c("MTurk", "Elite"), c(5,4)), 3),
         split_quantitle = rep(c(50, 25, 75), each = 9),
         entry = paste0(format_num(estimate, 2), " ", add_parens(std.error, 2)))


agree_df <-
  fit_df %>%
  dcast(term ~ split_quantitle + study, value.var = "entry") %>%
  mutate(term = c("Amtrak", "Climate", "Flat Tax", "Veterans", "Wall Street"))


summary_row <-
  fit_df %>%
  group_by(study, split_quantitle) %>%
  summarise(summary_est = weighted.mean(estimate, 1 / std.error ^ 2),
            summary_se = sqrt(1 / sum(1 / std.error ^ 2))) %>%
  mutate(entry = paste0(format_num(summary_est, 2), " ", add_parens(summary_se, 2)),
         term = "Precision weighted average") %>%
  dcast(term ~ split_quantitle + study, value.var = "entry")

print(xtable(agree_df), include.rownames = FALSE, include.colnames = FALSE, hline.after = c(), only.contents = TRUE)
print(xtable(summary_row), include.rownames = FALSE, include.colnames = FALSE, hline.after = c(), only.contents = TRUE)

